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A Generalized Labeled Multi-Bernoulli Filter 
Implementation using Gibbs Sampling 

Hung Gia Hoang, Ba-Tuong Vo, Ba-Ngu Vo 


Abstract —This paper proposes an efficient implementation of 
the generalized labeled multi-Bernoulli (GLMB) filter by com¬ 
bining the prediction and update into a single step. In contrast 
to the original approach which involves separate truncations in 
the prediction and update steps, the proposed implementation 
requires only one single truncation for each iteration, which 
can be performed using a standard ranked optimal assignment 
algorithm. Furthermore, we propose a new truncation technique 
based on Markov Chain Monte Carlo methods such as Gibbs 
sampling, which drastically reduces the complexity of the filter. 
The superior performance of the proposed approach is demon¬ 
strated through extensive numerical studies. 

Index Terms —Random finite sets, delta generalized labeled 
multi-Bernoulli filter 


I. Introduction 

Multi-object tracking refers to the problem of jointly es¬ 
timating the number of objects and their trajectories from 
sensor data. Driven by aerospace applications in the 1960’s, 
today multi-object tracking lies at the heart of a diverse range 
of application areas, see for example the texts d-ia. The 
most popular approaches to multi-object tracking are the joint 
probabilistic data association filter la, multiple hypothesis 
tracking 0, and more recently, random finite set (RFS) 0, 
0 . 

The RFS approach has attracted significant attention as 
a general systematic treatment of multi-object systems and 
provides the foundation for the development of novel filters 
such as the Probability Hypothesis Density (PHD) filter 1^- 
0, Cardinalized PHD (CPHD) filter 0, HO), and multi- 
Bernoulli filters 0, im, ca. While these filters were not 
designed to estimate the trajectories of objects, they have 
been successfully deployed in many applications including 
radar/sonar ca, m, Ea, computer vision Ca-llia, cell 
biology ca, autonomous vehicle Ea-Ea, automotive safety 
Ea, Ea, sensor scheduling Ebl - Ea and sensor network 

oni-oa- 

The introduction of the generalized labeled multi-Bernoulli 
(GLMB) RFS in has led to the development of the first 
tractable RFS-based multi-object tracker - the (5-GLMB filter. 
The 5-GLMB filter is attractive in that it exploits the conjugacy 
of the GLMB family to propagate forward in time the (labeled) 
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multi-object filtering density exactly ES- Each iteration of 
this filter involves an update operation and a prediction oper¬ 
ation, both of which result in weighted sums of multi-target 
exponentials with intractably large number of terms. The first 
implementation of the ^-GLMB filter truncate these sums by 
using the iL-shortest path and ranked assignment algorithms, 
respectively, in the prediction and update to determine the most 
significant components 

While the original two-staged implementation is intuitive 
and highly parallelizable, it is structurally inefficient as it 
requires many intermediate truncations of the J-GLMB den¬ 
sities. Specifically, in the update, truncation is performed by 
solving a ranked assignment problem for each predicted <5- 
GLMB component. Since truncation of the predicted J-GLMB 
sum is performed separately from the update, in general, a 
significant portion of the predicted components would gen¬ 
erate updated components with negligible weights. Hence, 
computations are wasted in solving a large number of ranked 
assignment problems, each of which has cubic complexity in 
the number of measurements. 

In this paper, we present a new implementation by formu¬ 
lating a joint prediction and update that eliminates inefficient 
truncation procedures in the original approach. The key inno¬ 
vation is the exploitation of the direct relationship between the 
components of the (5-GLMB filtering densities at consecutive 
iterations to circumvent solving a ranked assignment problem 
for each predicted component. In contrast to the original im¬ 
plementation, the proposed joint implementation only requires 
one truncation per component in the filtering density. 

Naturally, the joint prediction and update allows truncation 
of the (5-GLMB filtering density (without explicitly enumerat¬ 
ing all the components) using the ranked assignment algorithm 
03, ES, E3- More importantly, it admits a very efficient 
approximation of the (5-GLMB filtering density based on 
Markov Chain Monte Carlo methods. The key innovation is 
the use of Gibbs sampling to generate significant updated 5 - 
GLMB components, instead of deterministically generating 
them in order of non-increasing weights. The advantages 
of the proposed stochastic solution compared to the rank 
assignment algorithm are two-fold. First, it eliminates unnec¬ 
essary computations incurred by sorting the components, and 
reduces the complexity from cubic to linear in the number of 
measurements. Second, it automatically adjusts the number of 
significant components generated by exploiting the statistical 
characteristics of the component weights. 

The paper is organized as follows. Background on labeled 
RFS and the (5-GLMB filter is provided in section HI] Section 
Hin presents the joint prediction and update formulation and 
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the Gibbs sampler based implementation of the (5-GLMB filter. 
Numerical results are presented in Section |IV] and concluding 
remarks are given in Section [V] 


II. Background 


This section summarizes the labeled RFS and the GLMB 
filter implementation. We refer the reader to the original work 
lisa, da for detailed expositions. 

For the rest of the paper, single-object states are represented 
by lowercase letters, e.g. x, x while multi-object states are 
represented by uppercase letters, e.g. X, X, symbols for 
labeled states and their distributions are bolded to distinguish 
them from unlabeled ones, e.g. x, X, tt, etc, spaces are 
represented by blackboard bold e.g. X, Z, L, N, etc, and the 
class of finite subsets of a space X is denoted by J^(X). We use 
the standard inner product notation {f,g) = J f{x)g{x)dx, 
and the following multi-object exponential notation = 
Yixex ^i^)’ where is a real-valued function, with = 1 
by convention. We denote a generalization of the Kronecker 
delta that takes arbitrary arguments such as sets, vectors, etc, 
by 


6y{X) ^ 


1, if X = Y 
0, otherwise ’ 


and the inclusion function, a generalization of the indicator 
function, by 


Iy{X) ^ 


1, if X c y 

0, otherwise 


We also write 1^(3^) in place of li^({a:}) when X = {a;}. 


A. Labeled RFS 

A labeled RFS is simply a finite set-valued random variable 
where each single-object dynamical state is augmented with a 
unique label that can be stated concisely as follows 

Definition 1. A labeled RFS with state space X and (discrete) 
label space L is an RFS on XxL such that each realization 
has distinct labels. 


Let C : XxL —L be the projection C{{x,i)) = £, then 
a finite subset set X of XxL has distinct labels if and only 
if X and its labels £(X) = {i3(x) : x G X} have the same 
cardinality, i.e. b|x|(|-L(X)|) = 1. The function A(X) = 
(5|x|(|>L(X)|) is called the distinct label indicator. 

The set integral defined for any function f : X(XxL) —M 
is given by 

f f(X)bX==^i f f({xi, ...,xj)d(xi, ...,X,;). 
i=0 *■ 


where the integral of a function f XxL —>■ R is: 

J f{x)dx = Y^ J f{{xj))dx, 


eeh • 


The notion of labeled RFS enables the incorporation of 
individual object identity into multi-object system and the 
Bayes filter to be used as a tracker of these multi-object states. 


B. Bayes filter for labeled RFS 

Suppose that at time k, there are Nk target states 
Xfc.i,..., Xfc each taking values in the (labeled) state space 
X X L. In the random finite set formulation the set of targets 
is treated as the multi-object state 


Xfc = {Xfc4, . . . , Xfc_Ar^}. (I) 

Each state {xk,£) G Xfc either survives with probability 
Ps{xk,£) and evolves to a new state {xk+i,i) or dies with 
probability 1 — ps{xk,i). The dynamics of the survived 
targets are encapsulated in the multi-object transition density 
f/c-|-l|/c(X|Xfc). 

For a given multi-object state X^, each state {xk,£) G X^ 
at time k is either detected with probability PD{xk,i) and 
generates an observation z with likelihood g{z\xk, i) or missed 
with probability 1 — pjj{xk,£)- The multi-object observation 
at time k, Zk = {zk.i, •. •, Zk,Mk}^ the superposition of the 
observations from detected states and Poisson clutters with 
intensity k. Assuming that, conditional on Xfc, detections are 
independent, and that clutter is independent of the detections 
and is distributed as a Poisson RFS, the multi-object likelihood 
is given by ||33l, 041 

5(Zfc|Xfc) = ^ (2) 


where 0: L ^ {0,1,..., \Zk\} is a function such that 6{i) = 
9{i') > 0 implies i = i', and 


i; 0) = 


K-(ze{t)) ’ ^ I 

1 -Pd{x,£), if 6>(f) = 0 


(3) 


Remark. 6 is called an association map since it provides the 
mapping between tracks and observations, i.e. which track 
generates which observation, with undetected tracks assigned 
to 0. The condition 9{i) = 9(i') > 0 implies i = i' ensures 
that a track can generate at most one measurement at a point 
of time. 

The image of a set / C L through the map 9 is denoted by 
9{I), i.e. 

9{I) = {9{i) ■.£€!}, 


while the notation 0/ is used to denote the collection of all 
eligible association maps on domain I, i.e. 

Qi = {9-.9-\{0,l,...,\Zk\})=l} 

If the clutter is distributed as an iid cluster RFS, i.e. 
TTciK) = p{\K\)\K\\[pc{-)\^, the multi-object likelihood is 

g{Zk\Xk)=Y,^c{Zk\Zc) Y. 9M{Zc\Xk-,9c), (4) 

where 9q is the association map from £(Xfc) to {0,..., \Zq\} 
and 

Given a multi-object system as described above, the ob¬ 
jective is to find the multi-object filtering density, denoted 
by 7rfc+i(X|yfc+i fl which captures all information on the 


* For convenience, we drop the dependence on past measurements upto time 
k 
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number of targets and individual target states at time fc +1. In 
multi-object Baysian filtering, the multi-object filtering density 
is computed recursively in time according to the following 
prediction and update, commonly referred to as multi-object 
Bayes recursion m 


— J ffc+i|fe(X|Xfc)7rfc(Xfc|Zfe)5K/j, 

g{Zk+l\^''^k+l\k{^\Zk} 


7rfe+i(X|Zfc+i) = 


/5(Zfe+i|X)^fc+i|fc(X|Zfe)5X- 


(5) 


( 6 ) 


Note, however, that the Bayes filter is intractable since the 
set integrals in (HI)-® have no analytic solution in general. 


/(!) = 

rtf(i) /(i) T 

5 • • • 

/( 2 ) = 

r/( 2 ) ^( 2 ) 1 

1^1 , ■ • • 7{2)|/ 


/('*) = 

tll(k) M 1 


^( 1 ) 

^( 2 ) 









pd)(.,4i)) 



p('*)(-,4'‘’) 



Fig. 1. An enumeration of a 5-GLMB parameter set with each component 
indexed by an integer h. The hypothesis for component his ) while 

its weight and associated track densities are and {■,£), I G . 


C. Delta generalized labeled multi-Bernoulli RFS 

The i5-GLMB RFS, a special class of labeled RFS, provides 
an exact solution to ®-®. This is because the 5-GLMB 
RFS is closed under the multi-object Chapman-Kolmogorov 
equation with respect to the multi-object transition kernel 
and is conjugate with respect to the multi-object likelihood 
function ||3^ . 


Definition 2. A 5-GLMB RFS is a labeled RFS with state 
space X and (discrete) label space L, distributed according to 

X 


,(€) 


(7) 


^(X) = A(X) ^ 

(/. 5 )e.F(L)xH 

where S is a discrete space while and satisfy 

^ = 1 , ( 8 ) 


(/.«)GJr-(L)xE 




(9) 


Remark. The (5-GLMB density is essentially a mixture of 
multi-object exponentials, in which each components is iden¬ 
tified by a pair (/,^). Each I G -FjL) is a set of tracks 
labels while ^ € S represents a history of association maps 
^ = {di,... ,9k). The pair (/,^) can be interpreted as the 
hypothesis that the set of tracks / has a history of ^ association 
maps and corresponding kinematic state densities The 
weight uj^^’^^6i{£{lC)), therefore, can be considered as the 
probability of the hypothesis (/,^). 

Denote the collection of all label sets with n unique 
elements by J^„(L), the cardinality distribution of a 5-GLMB 
RFS is given by 


p{n) = ^ ( 10 ) 

(f.?)e3Pn(IL)xH 

A (5-GLMB is completely characterized by the set of param¬ 
eters : (1,5) € Jj^L)x5}. For implementation 

it is convenient to consider the set of parameters as an enu¬ 
meration of all 5-GLMB components (with positive weight) 
together with their associated weights and track densities 

as shown in Figure [T] where 

^(h) A andpW 

Given a 5-GLMB initial density, all subsequent multi-object 
densities are (5-GLMBs and can be computed exactly by a 
tractable filter called the ^-GLMB filter. 


D. Delta generalized labeled multi-Bernoulli filter 

The (5-GLMB filter recursively propagates a (5-GLMB den¬ 
sity forward in time via the Bayes recursion equations ® and 
®. Closed form solutions to the prediction and update of the 
()-GLMB filter are given in the following propositions 1^ . 

Proposition 1. If the multi-target posterior at time k — 1 is a 
5-GLMB of the form ®, i.e. 


7rfe_i(X|Xfc_i) = 


(4-1.&-i) 


AX) E 

( 4 -i.&-i)£. 4 (L)x: 


SiU^iX)) pt-A■\Zk.^) 


X 


( 11 ) 


and the birth density ts is defined on J^(X x B) according to 


fB(Y) = A(Y)wB(/:(Y))[pB(-)r, (12) 


then the multi-target prediction density to the next time is a 
5-GLMB given by 


'Tfc|fc-i(X|^fe-i) — A(X) X 


( 4 , &- 1 ) £ 4 '(LU B) X! 




iX 


,(13) 


where 


n L)u;B(/fc n B) 




(L) = 


Bs 


Xt 

4-i3i 




4-1-i 


Pklt-li^’ f) + iKitjPBix, £) 


pf ^\x,£) = 


Psi-J)fk\k-i{x\-,i),Pk!li\-,£) 


(14) 

(15) 

(16) 

(17) 

(18) 


Proposition 2. Given the prediction density in (113b . the multi¬ 
target posterior is a 5-GLMB given by 


7Tk{X\Zk) = A(X) 


E 






(Zk) X 


( 4 .&-i.Sfc) £ 4 (LUB) X H X 0 J-J, 

^4(-C(X)) 


pf-^’^’‘\-\Zk) 


(19) 






















4 


where 




{Zk) 


(&-i,4)/ ^'fc) 

P. 


{h,ik-l) 

‘^k\k-l 

J&-i) 




( 20 ) 

( 21 ) 

( 22 ) 


The propagations of (5-GLMB components through predic¬ 
tion and update are illustrated in Fig. |2] and Fig. [2 respectively. 
It is clear that the the number of components grows expo¬ 
nentially with time. Specifically, a component in the filtering 
density at time k — 1 generates a large number of predicted 
components, of which each one in turn produces a new set of 
multiple ()-GLMB components in the filtering density at time 
k. Hence, it is necessary to reduce the number of b-GLMB 
components in both prediction and update densities at every 
time step. 



7-(2) _ 


ri-^) _ 



(4'’. ■■■, 44 ^ 




Ak) 1 


f7i 

t(2) 

^fc-1 


Ah) 

^k-1 


^ jH 

1 



UJ 

{h) 

fc|fc-l 


„(1) f ^(1)\ 

„(2) , /2). 
Pklk-lV,H > 


(h) 

Pk\k- 

A;A) 


y(l) f /I) 

„{2) / .(2) ^ 


Ah) 

Pk\k-1 

(■ ^ 


/l\ 




II 

rW _ 


Ah) _ 







Ak) 1 


M) 

sfe 

.(h,2) 

^k 



{h,3) 

k 



^k 


, ,{h,j) 

‘^k 





pA\;A) 


(^.1)/ p{h) s 

Pk 

Ah) , 

Pk LVPV 


(hj) 

Pk 

(■ I 



j(i) _ 

11. . |^(i_)^|/ 

,(2) 

-'fc-l ^ 

11. . 


Ah) 

^k-1 — 

p{h) X 

’icA 


.(1) 

sfc-i 

^(2) 

Sfe-1 


Ah) 

Sfc-l 


“f-i 



47 


p7i(-.47 

pA{;A) 


pfA,A) 


„(i) (. /(i) 1 

„(2) (. A2) n 


pAAAa 


/l\ 

- 



j{h,l) _ 

AK2) _ 


j{h,j) _ 



!pi.h,\) p{h,l) \ 

Pi 1) 11 

11’ ’ 



p(h,j) \ 

i)|l 


Ah) 

^k-1 

Ah) 

^k-1 


Ah) 

^k-1 



, ,lk,2) 

“ilfc-l 


, ihj) 

“lc|i!-l 


pi^^-i(-.4"’^’) 

Ah) / p{h,2)\ 
Pk\k-A Ai } 


J.h) / p{h,j), 

Pfc|/t-l'. .‘1 ) 


Ah) ( piKi) s 
Pk\k-1\ 

(h) / pih,2) . 


Ah) ( p{h,3) \ 

Pk\k-1\ A^j(h,j)^) 



Fig. 2. The 5-GLMB prediction 1^ : component h in the prior generates 
a large set of predicted components with C U B, i.e. j = 1,..., 


and 

k \ k—1 k I fc—1 


.{h,3) ^(h)\ 


The simplest way to truncate a i5-GLMB density is dis¬ 
carding components with smallest weights. The following 
proposition asserts that this strategy minimizes the Li-distance 
between the true density and the truncated one ll34l 


Fig. 3. The 5-GLMB update (34): component h in the 
generates a (large) set of update components with ^ 


■ . ^ (h,j) A Vk ^ 1^ 

weights ul ,3 = 


{h) I. 


predicted density 
and 


be an unnormalized 5-GLMB density, /f T C H then 


life — frill 


fn 

fr 

llfnlli 

fr 1 


= E - 


{1,0 


< 2 


(/.5)gH-T 

llfnlli - llfrlli 


llfclli 


III. Joint prediction and update eor the (5-GLMB 

FILTER 

In this section, we briefly review the original implementa¬ 
tion of the (5-GLMB filter in subsection IIII-AI and propose a 
new implementation strategy with joint prediction and update 
in subsection IIII-BI 


A. The original implementation 

The first implementation of the b-GLMB filter, detailed in 
IIM), recursively calculates the filtering density by sequentially 
computing the predicted and update densities at each iteration 
based on Proposition [T] and Proposition |2] Since direct imple¬ 
mentation of equations (fTSl l and (fT9T l is difficult due to the 
sum over supersets in (flSl l. the predicted and update densities 
are rewritten as (|23l) and respectively, with 


,(4- 




h-i-L 


(&-i) 

Bs 




Proposition 3. Let ||f||j^ = f |f(X)| SX. denote the Li-norm 
off : J^(XxL) —>■ K., and for a given IH CJ^(L) x S let 


fH(X) = A(X) ^ (u(^’«)(5/(£(X)) 


(F.f)eH 


In the prediction stage 123, each component {Ik-i, Cfc-i) 
with weight 4-7 generates a set of prediction compo¬ 
nents (L U J, ^k-i ) with weight 


(LUJ,5fc_i) 

^k\k-l 


= W, 


(7fc-i,Cfc-i), (7fe_i,Cfc-i), 


‘k-l 


UJc 


\L)oob[J), 
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7rfc(X|Zfe) = ■ 


h-i,Ck-hL,J,9k 


' i^k) 


LUJ 


&uj(£(X))[pf-"4-|^fc) 


E 


h-i-,^k-i:L,J,9k 


„w 


.liUj 


(23) 


(24) 


where L and J represent two disjoint label sets for survival and 
birth tracks, respectively. Since the weight of the prediction 
component can be factorized into two factors, 
and ujb{J), which depend on two mutually exclusive sets; 
truncating the predicted density is performed by solving two 
separate AT-shortest path problems for each set of tracks. 
This is because running only one instance of the iT-shortest 
path based on the augmented set of existing and birth tracks 
generally favours the selection of survival tracks over new 
births and typically results in poor track initiation. 

In the update stage (l24li . each prediction component {L U 
J,^k-i) generates a (large) set of update components (L U 
J, 0fc)). These update components are truncated without 
having to exhaustively compute all the components by solving 
a ranked assignment problem. 

Although the original two-staged implementation is intuitive 
and highly parallelizable, it is has several drawbacks. First, 
since truncation of the predicted 5-GLMB density is performed 
separately from the update based purely on a priori knowledge 
(e.g. survival and birth probabilities), in general, a significant 
portion of the predicted components would generate updated 
components with negligible weights. Hence, computations are 
wasted in solving a large number of ranked assignment prob¬ 
lems, each of which has cubic complexity in the number of 
measurements. Second, it would be very difficult to determine 
the final approximation error of the truncated filtering density 
as the implementation involves least three separate truncating 
processes: one for existing tracks, one for birth tracks, and one 
for predicted tracks. 

In the following subsections, we will introduce the joint 
prediction and update as a better alternative to the original two- 
staged approach. The joint strategy eliminates the need for sep¬ 
arate prediction truncating procedures, thus involves only one 
truncation per iteration. Consequently, the new implementation 
yields considerable computational savings while preserving the 
filtering performance as well as the parallelizability of the 
original implementation. 


B. The joint prediction and update implementation 

Instead of computing the filtering density in two steps, 
the new strategy aims to generate the components of the 
filtering density in one combined step by formulating a direct 
relationship between the component of the current filtering 
density with those of the previous density. Specifically, we will 
derive a new formulation for 7rfc(X|Zfc) that does not involve 
prediction induced variables L and J. This can be done via 


an extended association map, denoted by 9, and defined as 
follows 


Definition 3. The extended association map is a function 9 : 
L U B —> {0, 1,... ,\Zk\,\Zk\ + \j such that 9{i) = 9{i') for 
Q<9{i) < |Zfc|-|-l implies i = i'. 


Remark. The new map, in essence, only extends the original 
map to include a new association, |Za;|-|- 1. In particular, 9 is 
identical to 9 except for non-survival and unconfirmed birth 
tracks, i.e. 


m 


9{e) We LUJ, 

\Zk\ + l V£ e (4_1 - L) U (B - J). 


The image of a set / C L U B through the extended map 
9 and the collection of all eligible extended association maps 
on domain I are denoted by 9{I) and 0/, respectively. 

Based on the notion of extended association map, the fol¬ 
lowing proposition establishes the direct relationship between 
two consecutive filtering densities at time k and k — 1. For sim¬ 
plicity, we assume that target births are modeled by (labeled) 
multi-Bernoulli RFS’s, i.e. ujb{J) = [1 ~ with 

r{t} denotes the existence probability of track f. 


Proposition 4. If the multi-target posterior at time k — \ is 
a 5-GLMB of the form (E]) and the set of targets bom at the 
next time is distributed as a labeled multi-Bernoulli RFS, then 
the multi-target posterior at the next time is a 5-GLMB given 
by 


7r,(X|Z,)=A(X) ^ 




i{o.....iz.|}(0fc(/:(x)))U«'=-"4l^fc) 


-iX 


(26) 


where Vf: 9k{f) < \Zk\ + l and 


, ,(4-l,&-l,4) , (Ik-1,ik-l) 

* W-1 




-i/fc-iUB 




1 - W G /fc_i :9k(i) = \Zk\ + l, 

l-r(f) VfGB:0fc(£) = |Zfe| + l, 

ri()ilzk^() V£GB:0fc(f)<|Zfc| + l. 


(27) 


(28) 


We now proceed to detail an efficient implementation of 
the i5-GLMB filter based on the result in Proposition 01 Let 
the sets of existing tracks, birth tracks, and measurements be 
enumerated by = {£i,.. .,£n}, B = {In+i, ■ ■ -Ap}, and 
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Zk = {^1, ■ • ■, zm}, respectively. Given the /i-th component 
(^k-h ^k-i)’ objective is to find extended associations 

J produce highest update weights. 

Denote by Sj a P x [M + 2P) matrix whose entries are 
either 1 or 0 such that the sum of each row is exactly 1 
while the sum of each column is at most 1. The matrix Sj is 
called an assignment matrix since it represents a valid extended 
association map G Hence, finding the desirable 

extended associations is equivalently translated to finding 
the corresponding assignment matrices. The simplest way to 
determine these matrices without exhaustingly computing all 
of the update weights is to assign an appropriate cost for 
each assignment matrix and then rank these matrices in non¬ 
decreasing order of their costs via Murty’s algorithm. 

Let be a matrix whose (n, m) entry, with 1 < n < P 
and 1 < m < M -f 2P, is defined as follows. 


’(4) 




i^n) 

.0 


m < M, 
m = M -\- n, 
m = M + P + n, 

otherwise. 


(29) 


As depicted in Fig. |4] the matrix F^^ contains all possible 

values of 72^*^ {£) that any valid extended associations 

can generate. Intuitively speaking, each entry of F^^ is an 
indicator of how likely an extended association is assigned 
to a track; hence, the matrix = — log ^F^^ ^ is called 
the cost matrix and cost function for a particular assignment 
matrix Sj is given by 

f4S,) = -tr[sfCp!), (30) 

It is straightforward to show that the (unnormalized) weight of 
the -generated component in the filtering density at time 
k is 

ujk (xexp[-fc[Sj)]. (31) 

Based on (l30l l and OTl i. an implementation of the 5-GLMB 
filter with joint prediction and update is given in Algorithm [T] 
where the subroutine ranked_assignment uses Murty’s algo¬ 
rithm ll^ to generate a sequence of assignment matrices 
that yield lowest costs (or equivalently, highest update weights) 
without exhaustively navigating the whole assignment space. 


Similar to the original implementation, the joint prediction 
and update also operates independently on each components 
in the filtering density, thereby is highly parallelizable. 


C. Stochastic simulation approach to extended data associa¬ 
tions 

In multi-target tracking, truncating procedures based on the 
original Murty’s algorithm with complexity o(t\Z\‘^, where 
T is the number of assignments and \Z\ is the number of 
the measurements, have been proposed in lEsi, mu. Ho). 
More efficient algorithms with 0\T\Z\j complexity have 


Algorithm 1 i5-GLMB joint prediction and update 


Inputs: 

Outputs; {^ j 


H 

Zk 


1 : for t— 1, do 

ih) 


(Li)=(i.i) 


compute F^^^ according to (??) 

cf> ^ - i4(r|) 

■G- ranked_assignment(C2^\ 
for j ^ 1, rW do 


compute according to (EB 

compute according to (l22l i 


end for 
end for 


10: normalize 


weights 




been proposed in m, Il36l, EH, with the latter showing 
better efficiency for large \Z\. The main drawback of these 
approaches is that a significant amount of computation is 
used to sort the data associations in a particular order despite 
the fact that order is effectively discarded after the update. 
Furthermore, the number of desired components must be 
predetermined, generally by a large enough number to capture 
all important associations. If the number of significant compo¬ 
nents in the filtering density is much smaller than the chosen 
threshold, many insignificant components are generated that 
waste a lot of computation at the next iteration. Conversely 
if the number of significant component exceeds the chosen 
threshold, the filtering performance will likely degrades in 
subsequent iterations. Nonetheless a deterministic polynomial 
time solution is thus appealing in the sense that convergence 
and reproducibility is guaranteed without having to enumerate 
all possible solutions. 

In this subsection, we propose an alternative to the ranked 
assignment based solution. Our proposed solution is based on 
stochastic simulation or Markov Chain Monte Carlo methods 
via a Gibbs sampler which directly address the above men¬ 
tioned drawbacks. Conceptually, instead of ranking extended 
associations in non-increasing order of their weights, each 
extended association is treated as a realization of a (discrete) 
random variable, where the probability of each extended 
association is proportional to the weight of its associated i5- 
GLMB component in the next filtering density. Candidate 
extended associations are then generated by sampling from this 
discrete distribution. Extended associations with high weights 
are chosen more often than those with low weights in a 
statistically consistent manner. Consequently these samples 
of extended associations are more statistically diverse than 
those from obtained from a deterministic approach such as 
the ranked optimal assignment. 

Using the same enumeration for tracks and measurement 
as in the previous section, a valid extended association map 9 
for each component can be represented as a vector 

r.. .. tT 


9 = 


9{ii),...,9iip) 


G {0,1,..., M-b 1}^. The key idea 
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Fig. 4. The matrix F^ whose entry represents the likelihood that = m for any valid extended association The cost 

■^fc—1 

matrix is foiTned by taking negative logarithm on , i.e. = — log ^. 


of the stochastic based approach is that 9 can be considered 
as realizations of a random variable in the space 0 ,(»>), with 

r fc_l UB 

the following distribution 


7r(0) oc 






(i(^) 


p(h) 




(32) 


where 




(»)=n [ 


1 -1 


/r(^) 

fc _ 1 tSfc — 1 5^ j 


2=1 



(33) 

) r (ft (')') lUB 

(34) 


with M - n {0(4),..., 0>p)} and 

is given in (l28l l. Thus the probability of a valid extended 
association is proportional to the weight of the corresponding 
(5-GLMB component in the next filtering density while zero 
probability is allocated to extended associations which do not 
satisfy the constraint that each measurement is assigned to at 
most one track. 

However, sampling directly from the distribution (l32T i is 
very difficult since we cannot exhaustively compute all of 

the values of .A common solution to this kind 

of problem is to use Markov Chain Monte Carlo (MCMC) 
methods such as the Gibbs sampler to obtain samples from 

(|32] | without having to directly compute . The 

Gibbs sampler is a very efficient method to sample a difficult 
distribution if its conditional marginals can be computed in a 
simple closed form Il42l . Il43l with proven convergence under 
generally standard assumptions ia, ia. 

The main theoretical contribution in this section is stated in 
the following proposition, which allows conditional marginals 
to be computed via the entries of the matrix 


Proposition 5. Denote by 9n the n-th element of 6, i.e. 
Sn = 0(^ri), ond Ofi all the other elements except On, i.e. 
0n = [01:n—1, 0n+l:p] . Then, the conditional marginal 7r(0„|0fl) 
is given by 


1- 

r(0„|0fi) oc 1^: 


1-1 






(35) 


i=l 

Remark. Propostion [3 simply states that the conditional prob¬ 
ability 7r(0„|0fi) are zero when 0„ is inconsistent with 9m, thus 


ensuring that each measurement can be assigned to at most on 
track. The conditional probabilities are otherwise proportional 

to itf’iQ- 

With these conditionals the pseudo code for (5-GLMB filter¬ 
ing with the Gibbs sampler is given in Algorithm where the 
Gibbs sampler is used directly in place of the ranked optimal 
assignment. In order to produce one sample. Algorithm |2] 


Algorithm 2 Gibbs sampling of (5-GLMB assignments 
Inputs: il, U B, 

j'Qt) 


Outputs: 






2: compute according to 
3: for j ^ 1,T('*) do 


Initialize [si,..., sp]^ 

for f 1, do 

for z <— 1, P do 

compute 7r(si|si:i-i, Si+i:p) according to {33 
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13: 
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ih,t). 


7r(si|si:i-l, Si-|-l:p) 


15: end for 


end for 
end for 
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fh,j) Ah,j)\ 


requires a Gibbs sequence of length starting from an arbi¬ 


trary initialization 


fhfi). 


’^(4),...,0f%p) 

Alternatively, we can sample a long Gibbs sequence of length 
B(t‘) X + 1 and then extract every P^^^-th sample Bbll . 
The length of the Gibbs sequence, roughly speaking, depends 
on the convergence rate of the Gibbs sampler and the distance 
from the initial point to the true sample space. If we start with 
a good initialization right in the true sample space, we can 
use all the samples from the Gibbs sequence 1471 . In practice, 
one example of good initialization that allows us to use all of 


, to be generated. 








































Fig. 5. Multiple trajectories in the xy plane. Start/Stop positions for each 
track are shown with o/A. 


the samples from the resulting Gibbs sequence is the optimal 
assignment, which can be obtained via either Munkres ll48l 
or Jonker-Volgenant algorithm 1491 . Otherwise, we can start 
with all zeros assignment (i.e. all tracks are misdetected) that 
is also valid sample and requires no additional computation. 

In terms of computational complexity, sampling from a 
discrete distribution is linear with the weight’s length l50l . 
therefore the total complexity of the Gibbs sampling procedure 
presented in Algorithm |2] is P{M + 2P)). In 

comparison, the fastest ranked optimal assignment algorithm is 
0{T^^\M-\-2PY) 1^ . IJTI . For general multi-target tracking 
problems in practice, we usually have P <C (M + 2P), 
thus the Gibbs sampling algorithm will generally be much 
faster than the ranked assignment given the same 

IV. Simulation 


transition density /fe|fe_i(a;fc|a;fc-i) = M{xk\FkXk-i,Qk) and 
likelihood giJyZk\xk) = Af{zk] HkXk, Rk) with parameters 


Fk 


I2 A/2 
O 2 h 


Hk 


h O 2 


Qk — O'jy 
Rk = 


^ T 

—h 

^ T 
— ^2 

T 

— ^2 

A2/2 


where /„ and 0„ denote the n x n identity and zero 
matrices respectively, A = Is is the sampling period, 
cTi, = 5m/s‘^ and = lOm are the standard deviations 
of the process noise and measurement noise. The survival 
probability is ps,k = 0.99 and the birth model is a La- 

(ll ('D ^ 

beled Multi-Bernoulli RFS with parameters ttb = f’B^PBi=i 
where r® = 0.04 and p®(a;) = A/’(a;; m®, Pg) with 

TO® = [0,0,100,0]"^, = [-100,0,-100,0]'^,TO*| = 

[100, 0, -100,0]^, Pb = diag([10,10,10,10]^)^ . The detec¬ 
tion probability is pD,k = 0.88 and clutter follows a Poisson 
RFS with an average intensity of Ac = 6.6 x lCr®m“^ giving 
an average of 66 false alarms per scan. 

First, we compare the performance of the traditional sepa¬ 
rated and the proposed joint prediction and update approaches. 
For a fair comparison, both approaches are capped to the same 
maximum components. Results are shown over 100 Monte 
Carlo trials. Figures |6] shows the mean and standard deviation 
of the estimated cardinality versus time. Figures [T] and [8] show 
the OSPA distance iM] and its localization and cardinality 
components for c = 100m and p = 1. 


Cardinality statistics for joint prediction and update 



In this section we first compare the performance of the joint 
prediction and update approach with its traditional separated 
counterpart, both employ the ranked assignment algorithm for 
fair comparison. Then, we illustrate the superior performance 
of the Gibbs sampler based truncation to the conventional 
ranked assignment via a difficult tracking scenario with low 
detection probability and very high clutter rate. 

The first numerical example is based on a scenario adapted 
from II 34 II in which a varying number targets travel in straight 
paths and with different but constant velocities on the two 
dimensional region [—1000, lOOOjm x [—1000, lOOOjm. The 
duration of the scenario is K = lOOs. There is a crossing 
of 3 targets at the origin at time k = 20, and a crossing of 
two pairs of targets at position (±300,0) at time k = 40. The 
region and tracks are shown in Figure |5] 

The kinematic target state is a vector of planar position 
and velocity Xk = [Px,k,Py,k,Px,k,Py,k]'^ ■ Measurements are 
noisy vectors of planar position only Zk = [Zx,k, Zy^k]"’" ■ The 
single-target state space model is linear Gaussian according to 


Cardinality statistics for separated prediction and update 



Fig. 6. Cardinality statistics with traditional ranked assignment truncation. 

It can be seen that both approaches estimate the cardinal¬ 
ity equally well. Similarly, in terms of OSPA distance, the 
performance of the two approach is virtually the same. 

Second, we demonstrate the fast implementation via the 
Gibbs sampler. In this example, we keep all parameters the 
same as in the previous example except that the clutter rate 
is now increased to average 100 false alarms per scan. The 
performance of the Gibbs sampler implementation is compared 
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Fig. 7. OSPA distance with traditional ranked assignment truncation. 


Localization error 



Cardinality error 



Fig. 8. OSPA components with traditional ranked assignment truncation. 



Fig. 9. OSPA distance compaiison between the Gibbs sampler implementa¬ 
tion versus the ranked assignment implementation (c = 100m, p = 1). 


with that of a ranked assignment based implementation with 
the same maximum number of posterior hypotheses. The 
average OSPA distances over 100 Monte Carlo trials are 
presented in Fig. |9] 

It is obvious that the Gibbs sampler has a better OSPA 
from around time k = 75 onward. The reason is in difficult 



Fig. 10. Cardinality statistics for Gibbs sampling implementation. 



Fig. 11. Cardinality statistics for ranked assignment implementation. 

scenario (e.g. high clutter rate, low detection probability), if 
the number of existing targets are high the Gibbs sampling 
technique is expected to pick up the new born target better 
than the ranked assignment algorithm given the same number 
of samples/hypotheses due to its randomized behaviour. This 
is clearly illustrated in the cardinality statistics for both ap¬ 
proaches in Fig. [To] and Fig. [TT] As expected, however, the 
joint approach averaged run time is significantly lower that 
that of the original approach. Reductions in execution time of 
1 to 2 orders of magnitude are typical. 

V. Conclusions 

In this paper we propose a new implementation scheme for 
the b-GLMB filter that allows joint prediction and update. 
In contrast to the conventional two-staged implementation, 
the joint approach use a posteriori information to construct 
cost matrices for every individual track, thereby requires only 
one truncation in each iteration due to the elimination of 
inefficient intermediate steps. More importantly, this joint 
strategy provides the platform for the development of an 
accelerated randomized truncation procedure that achieves 
superior performance as compared to that of its traditional 
deterministic counterpart. The proposed method is also ap¬ 
plicable to approximations of the i5-GLMB filter such as the 
labeled multi-Bernoulli (LMB) filter 1521 . 
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